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PDF approach for turbulent scalar 
field: some recent developments 


By Feng Gao 


1. Motivation and objectives 

The probability density function (PDF) method has been proven a very useful 
approach in turbulence research. It has been particularly effective in simulating 
turbulent reacting flows and in studying some detailed statistical properties gener- 
ated by a turbulent field (see, e.g., Monin & Yaglom 1975, Bilger 1989, Pope 1985 
& 1990, Kraichnan 1990, Gao et al. 1992). 

There are, however, some important questions that have yet to be answered in 
PDF studies. Our efforts in the past year have been focused on two areas. First, a 
simple mixing model suitable for Monte Carlo simulations has been developed based 
on the mapping closure. Secondly, the mechanism of turbulent transport has been 
analyzed in order to understand the recently observed abnormal PDF’s of turbulent 
temperature fields generated by linear heat sources (Gollub et al. 1991, Jayesh & 
Warhaft 1991). — 

1.1 Needs for new mixing models in PDF approach 

It is well known that the PDF approach provides a closed form representation 
for the chemical reacting source terms, thus becoming a preferred choice among the 
available closure models for turbulent reacting flows (O’Brien 1980, Pope 1985). 
The general argument is that once the PDF of a reacting scalar, ip, is known, 
the mean reaction rate, which is the key quantity to be evaluated in reacting flow 
simulations, can then be readily calculated by 

(S) = J S(tp)P(ip)dip. 

Generally, two approaches have been commonly taken in applying the PDF method, 
namely, presumed PDF and full PDF methods. 

The presumed PDF method assumes a certain form for the PDF with some 
adjustable parameters, which are supplied by other conventional models (Bray 1980, 
Borghi 1988, Vervisch 1992). Restricted by the scope of the closure models used (for 
example k — e model), the parameters to be adjusted are generally limited to the low 
order statistics such as mean and variance, thus exposing an important shortcoming 
of the presumed PDF. In fact, one can construct different PDF’s P\ and P 2 that 
have the same mean and variance but totally different higher order statistics. If 
the reaction source term is highly nonlinear (which is true in most cases), the mean 
reaction rates obtained from Pi and P 2 may be very different, depending on which 
PDF form is chosen. 
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The natural way for obtaining the PDF is the full PDF method which simulates 
the PDF from its evolution equation. However, a major stumbling block in this 
approach has been the lack of a proper closure for the diffusion effect (O’Brien 
1980, Pope 1985). In order for a model to be accepted in practical simulations, 
it has to be physically reasonable and numerically easy to implement. Despite 
the theoretical success the mapping closure enjoys in treating the mixing effect in 
the PDF approach (Kraichnan 1990, Gao 1991, Pope 1991), it has been shown 
difficult and computationally intense to implement this closure in simulations (Gao 
&: O’Brien 1991, Valino & Gao 1992). The most commonly used model for diffusion 
effect in practical Monte-Carlo simulations remains the LMSE model (Pope 1992) 
which reads 

- <*)). (l) 

It is well known that this model does not relax a PDF. In fact, it can be easily 
shown that applying this model leads to 

= *■.«•). 

where 0' = 0 — (0). Therefore, the PDF so obtained can be very erroneous. This 
puts us in a rather awkward position. On one hand, we are attempting to use a 
highly sophisticated approach whose main promise is to provide accurate estimates 
for mean reaction rates. On the other hand, the PDF could be so contaminated 
that it does not reflect the true evolution of the fields being considered. 

It is, therefore, obvious that in order for the PDF method to live up to its 
promises, better mixing models which are easy to implement should be developed. 

1.2 PDF generated by non-uniform sources 

The experiments of Jayesh & Warhaft (1991) and Gollub et a/.(1991) indicate 
that the scalar PDF generated by a linear source term exhibits exponential tails. 
This result is rather surprising because it has generally been believed that the PDF 
so generated is a Gaussian distribution (Venkataramani & Chevray 1978, Tavoularis 
&; Corrsin 1981). This situation certainly demands a theoretical investigation. 

There are three basic processes that determine the distribution of a scalar PDF: 
the shape of the non-uniform source, the turbulent convection, and the molecular 
diffusion. A fluid particle leaving the source is convected by the turbulent velocity 
field to a certain observation point. Because of the chaotic nature of the velocity 
field, particles from different positions in the source all have certain possibilities 
of reaching the observation point, thus generating the fluctuations that reflect the 
characteristics of the source. In the absence of molecular diffusion, the PDF of the 
scalar is determined by the interaction of the source and the convection. 

There are a couple of reasons that justify the neglection of molecular diffusion 
in search for the mechanism of generating exponential tails. First, it is supported 
by theoretical arguments and numerical simulations that the molecular diffusion 
tends to relax a PDF to a Gaussian distribution. Although it has been shown 
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that the interaction between random turbulent advection and molecular diffusion 
distorts a Gaussian PDF to generate mild non-Gaussian tails (Gao et al. 1992), the 
clear exponential tails observed in these experiments cannot be explained within 
the frame of this interaction. Secondly, for high Reynolds (Peclet) number flows, 
the diffusive effect is very small in comparison with the turbulent transport (Taylor 
1935) which is responsible for bring around the fluctuations generated by the non- 
uniform source. 

Based on these arguments, the experiments mentioned earlier can be analyzed 
explicitly under some idealized conditions. This study suggests a mechanism which 
seems to provide an explanation for the observed tails. 


2. Accomplishments 

2.1 A mixing model for PDF simulations 
In dealing with turbulent reacting flow problems, it is generally accepted to sepa- 
rate the effects of mixing and reaction by time-splitting schemes. Since the reaction 
term is closed in PDF formulation, we will concentrate on proper modeling for the 
mixing effect. 

The mapping closure maps a known statistical field <f> (generally chosen as a 
multi- variate standard Gaussian field) to a surrogate field X(<f),t) whose statistics 
resemble those of the true scalar field ip (Chen et al 1989). Under homogeneous 
and non-reacting conditions, the mapping relation is governed by 


dX 

dt 


± dX 

u ( -^ + 


d 2 X 
d<j> 2 h 


( 2 ) 


where is determined by the time scale of turbulent mixing. It has been shown 
that this model provides an excellent representation for the mixing effect in the 
PDF approach (Gao 1991a, Pope 1991). 

In spite of its good physics, the mapping closure has posed great difficulties for 
numerical implementation, as pointed out earlier. The problem stems from the 
necessity that the fields be re-mapped at each time step, which is computationally 
intense. This problem worsens drastically as the dimension (number of the scalar 
quantities) involved increases. 

In the model we are proposing, the PDF is represented by a group of represen- 
tative “particles” which advance in time following certain laws, as is generally used 
in Monte- Carlo simulations. In case of reaction acting alone, for example, each 
particle is advanced by 

dip , cf , x 

ST = 


The task is, therefore, to develop similar models to describe the mixing effect. 

One way to generate such models is to use the mapping solution. It is well known 
that the general solution for (2) is (Gao 1991a, 1991b) 


384 


F. Gao 


where t = J* uj*dt and H„ are Hermite functions. The expansion coefficients a„ 
can be determined as 




1 ,d < n > X 


2 "/ 2 n! ' 


)«=o 


(4) 


Clearly, for reasonably well behaved mapping, a„ tends to zero rapidly. Therefore, 
we expect that a truncation of the right hand side of (3) at a relatively low order 


71 = 0 


X(<f>,t) = fl a nHn(^=)e- r 


(5) 


can approximate to a satisfactory degree of accuracy. 

A group of surrogate particles can be chosen according to (5) to represent the 
PDF, and each of these particles evolves according to 

11(4 + - W) = 0, (6) 

1=1 


where d/dr = d/(u*dt) and to* can be related to the scalar evolution time scale. In 
fact, if u? in 


dt 




can be provided by other models, such as the k — e model, it can be shown that 


UJ 


<± S^i 2 n n!p n _! 

2 n ■ 2”n!p„_i ’ 


(7) 


where fx a = (a Q +i/ai) 2 e 2ar . If the truncated form (5) is used, all fi Q where a >m 
should be set to zero. 

can be related to the certain order moments of the field considered. For 
example, if m = 3, it can be shown that 


v5(1+4/i? + 24^)V* 

12 -(- 144/i2 + 32/iJ H“ 864^2 

and 

_ F 4 (l + 4 + 24^2 ) 2 ~ 

48(1 + 48/i 2 + 2\&n\) 

Here F, axe moment coefficients as defined earlier and 


(8a) 

(8b) 


g(fii , fi 2 ) = 1 + 40 ii\ -f 336/i2 + b9h2^\^\ + 80/4 4- 17856/4- 

These equations can be solved iteratively, and our tests seem to suggest that /i a are 
very small and 

u?* « u?/2 


(9) 
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FIGURE 1. Evolution of scalar variance: DNS (dotted line); LMSE (solid line) 
and current model (dashed line), tu/l is the eddy-turn-over time. 
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FIGURE 2. Evolution of scalar flatness: DNS (dotted line); LMSE (solid line) and 
current model (dashed line), tu/l is the eddy- turn-over time. 

remains a good approximation. Taking the highly singular double-delta PDF as an 
example, it can be shown that /i 2or +i = 0 (because of symmetry) and fi 2 = 1/144, 
fi 4 = 1/25600, etc. 

It is noticed that if we choose m = 1, (6) recovers the equation in the LMSE 
model (1) and (7) shows that w* = lj/ 2. 

It should be pointed out that in the current model, we are only interested in the 
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Figure 3a. Initial scalar PDF. 
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FIGURE 3b. Scalar PDF at t\ujl = 0.178: DNS (dotted line); LMSE (solid line) 
and current model (dashed line). 

evolution of a group of surrogate particles whose statistics closely resemble those of 
the true turbulent field. These particles are generally not the fluid particles. 

Some tests have been conducted using the model with m = 3 and compared 
with the corresponding cases from DNS and the LMSE model. Figure 1 shows the 
evolution of scalar variance with a? matched from DNS data. It shows that (8) is 
indeed a good approximation. Figure 2 exhibits the evolution of the scalar flatness 
F\. While the LMSE model clearly does not relax the PDF, the current model 
catches on the trend of DNS. 
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FIGURE 3c. Scalar PDF at t 2 u/l = 0.685: DNS (dotted line); LMSE (solid line) 
and current model (dashed line). 



FIGURE 4. Sketch of the system considered. 


The plots of PDF evolution perhaps are more revealing. The initial PDF in the 
model simulations are generated according to that in DNS and normalized in the 
interval [0,1] (Figure 3(a)). The evolved PDF at two later moments are plotted in 
Figures 3(b) and 3(c). The improvement achieved by the current model is obvious. 

The accuracy of representing the PDF by a collection of surrogate particles de- 
pends on the shape of the PDF. Models with better accuracies can be developed 
along the same line by pushing m higher. However, it is expected that m = 3 should 
suffice for practical simulations. 
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2.2 PDF of temperature fields generated by a linear heat source 

As discussed earlier, the abnormal tails of a scalar PDF are mainly caused by 
the interaction between the non-uniform source and the random advecting velocity. 
Let r% and be the time scales of turbulent transport and molecular diffusion, 
respectively, it is well known that t^/t* ~ Re (Tennekes & Lumley 1973). For 
high Reynolds number flows, we simply assume that r<* — ► oo and r t — ► 0, namely, 
neglecting the molecular diffusion and assuming the velocity fluctuation is white in 
time. The consequences of these simplifications are explained elsewhere (Gao 1992). 

Hence, a particle detected at time t at the observation point (x, y ) can be traced 
back to an earlier time r when it was released from the source at (0, yo) and acquired 
the temperature T(0,y o ;r) (refer to Figure 4 for a sketch of the system being 
considered), i.e. 


x = U(t — r) + J u'dt 

It can be shown that (Gao 1992) 

x = U(t — r) + c 7 x yjt — rr\ 

where = 2u,L, and r; are standard independent Gaussian random variables. 
Here u* is the variance of velocity and Li the correlation length in the i direction. 
The PDF of temperature can be written as 

P(f;x,y,*) = (<5(T-r(0,y o ;r)) yOir < t , (12) 

where the average is taken over all possible particles released from the source prior 
to t. 

Consider the case of linear source, i.e, 

^(O^yo;^) = ay 0 . 


and y 


— y° = J u 


( 10 ) 


and y - yo = <r y Vt ~ rr 2 , (11) 


It can be shown that 


P(T;x,y,<) = F(T\x,y,t)G{T\x,y,t) 


where 


and 


2naa 


F ( T ) = nZZZ ex P[" _2 


r<trr\ 1 + 7 ) , n , 0 2 \ jt 

a{T) = w L ~ 7f>' ■ - exp( - £ /2a « • 


+ 41^7 


(13) 

(14) 

(15) 


Here, a — a x , ft = o x j(T y and 


J 2 = X 2 + 0 2 (— - y) 2 . 
a 


( 16 ) 



PDF approach for turbulent scalar field 


389 


It can be demonstrated that 


lim 
r— oo 


dlnG 

dT 


= 0, 


hence the main contribution to the tails comes from F(T) which decays exponen- 
tially as T — » oo. 

In linear theory, the fluctuation in x (i.e. r*) is neglected compare with U{t-r). 
This is equivalent to letting a — * 0 while keeping <r 9 = a/fi finite. In this limit, 
we have ft — ♦ 0. It can be easily shown that 


F(t) 


2naa. 


r U ,T , 21 
eXp ("W ( a ~ »> 1 


V 


and 


Hence, 



P(T) = 


Vu 


y/2xTTi 


X7r aa u 


r U ( T ^ 


which recovers what would be obtained using the traditional method. 

However, this approximation seems to be a rather awkward one because in tur- 
bulence, fi is generally not small. In fact, in a grid generated turbulence, u' ~ U 
and v 1 ~ 0 near the grid and u' ~ v' in the later stage. It doesn’t seem reasonable 
to neglect the fluctuation in x direction while keeping the terms of similar or lower 
order of magnitude in the y direction. 

For finite /?, we therefore should expect exponential tails under turbulent disper- 
sion at not too distant downstream locations. The exponent can be shown as 


A\yl* + £ 2 (- - y) 2 - *1 ~ ~^i T for T>1 - 

a 2 V a acr 2 

Another limit can be obtained if we look at far downstream locations and finite 
temperature range. In this limit, /? 2 ( J — y) 2 /x 2 <C 1, hence, 

-^1 - *1 ~ 

which also recovers the traditional result. This approximation is equally applicable 
if we have a strong mean scalar gradient. Therefore, if we move the probe to further 
downstream locations and/or increase the gradient of the scalar source, we will see 
more Gaussian-like distributions. 

Other quantities, such as mean and higher order moments, can be readily calcu- 
lated. The qualitative features predicted by this idealized model compare favorably 
with the available experiments. For detailed discussions, please refer to Gao (1992). 
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3. Future plans 

Given the recent exciting developments in large-eddy simulations of turbulent 
fields (see related reports in this volume), it seems only natural to extend this suc- 
cessful technique to the reacting flow problems. However, despite many attempts, 
the proper modeling of the reaction source terms remains an unsolved problem. 
Considering the wide range of different type of reactions encountered in engineering 
applications, the prospect of developing a suitable sub-grid-scale model for these 
source terms seems to be grim. The PDF approach can probably be usefully em- 
ployed to address this problem. 
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